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ABSTRACT 

At  present  the  preferred  method  for  fitting  a  general  curve  through  scat¬ 
tered  data  points  in  the  plane  is  orthogonal  distance  regression ,  i.e.,  by  min¬ 
imising  the  sum  of  squares  of  the  distances  from  each  data  point  to  its  nearest 
neighbour  on  the  curve.  While  generally  producing  good  fits,  in  theory  orthog¬ 
onal  distance  regression  can  be  both  biased  and  inconsistent:  in  practice  this 
manifest  itself  in  overfitting  of  convex  curves  or  underfitting  of  corners.  The 
paper  postulates  this  occurs  because  orthogonal  distance  regression  is  based 
on  an  incomplete  stochastic  model  of  the  problem.  It  therefore  presents  an 
extension  of  the  standard  model  that  takes  into  accounts  both  the  noisy  mea¬ 
surement  of  points  on  the  curve  and  their  underlying  distribution  along  the 
curve.  It  then  derives  the  likelihood  function  of  a  given  curve  being  observed 
under  this  model.  Although  this  cannot  be  evaluated  exactly  for  anything  other 
than  the  simplest  curves,  it  lends  itself  naturally  to  asymptotic  approximation. 
Orthogonal  distance  regression  corresponds  to  a  first  order  approximation  to 
the  maximum  likelihood  estimator  in  this  model:  the  paper  also  derives  a  sec¬ 
ond  order  approximation,  which  turns  out  to  be  a  simple  modification  of  the 
least  squares  penalty  that  includes  a  contribution  from  the  curvature  at  the 
closest  point.  Analytical  and  numerical  examples  are  presented  to  demonstrate 
the  improvement  achieved  using  the  higher  order  estimator. 
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Fitting  the  Most  Likely  Curve  through  Noisy  Data 


EXECUTIVE  SUMMARY 

Many  problems  in  computer  vision  and  elsewhere  centre  on  finding  some  member  of 
a  family  of  curves  that  best  fits  a  given  set  of  points  in  a  plane  (or,  more  generally,  on 
finding  a  member  of  an  m-dimensional  family  of  surfaces  that  best  fits  a  set  of  points  in 
n-dimensional  space).  At  present  the  preferred  method  for  fitting  a  general  curve  through 
scattered  data  points  in  the  plane  is  orthogonal  distance  regression,  i.e.,  by  minimising 
the  sum  of  squares  of  the  distances  from  each  data  point  to  its  nearest  neighbour  on  the 
curve.  Orthogonal  distance  regression  is  well  known  to  be  computationally  difficult,  and 
although  some  good  algorithms  are  now  available,  it  can  still  be  very  hard  to  numerically 
determine  the  optimal  fit.  This  paper  is  not  concerned  with  efficient  computation,  however. 
Rather  it  explores  the  more  fundamental  problem  of  the  flaws  in  the  stochastic  model 
on  which  orthogonal  distance  regression  is  based,  and  develops  an  improved  model  and 
associated  distance  measures  that  produce  better  fits.  While  generally  producing  good  fits, 
in  theory  orthogonal  distance  regression  can  be  both  biased  and  inconsistent:  in  practice 
this  manifest  itself  in  overfitting  of  convex  curves  or  underfitting  of  corners.  The  paper 
postulates  this  occurs  because  orthogonal  distance  regression  is  based  on  an  incomplete 
stochastic  model  of  the  problem. 

The  paper  stemmed  from  the  authors’  practical  experience  in  fitting  curves  to  ionogram 
data,  while  the  approach  taken  was  inspired  by  some  perceptive  remarks  by  Kanatani  on 
the  distinction  between  geometric  and  statistical  distance  measures  in  fitting  curves  to 
data.  It  was  found  that  curves  fitted  to  these  data  sets  by  orthogonal  distance  regression 
had  a  tendency  to  underfit  the  nose  of  the  trace  (i.e.,  not  to  extend  as  far  as  the  data). 
Theoretical  confirmation  of  this  bias  in  the  fitting  procedure  would  have  serious  implica¬ 
tions  for  ionogram  modelling  due  to  the  importance  of  accurately  locating  the  exact  extent 
of  the  nose  in  applications  of  interest. 

This  paper  therefore  presents  an  extension  of  the  standard  orthogonal  distance  regres¬ 
sion  model  that  takes  into  accounts  both  the  noisy  measurement  of  points  on  the  curve 
and  their  underlying  distribution  along  the  curve.  It  then  derives  the  likelihood  function 
of  a  given  curve  being  observed  under  this  model.  Although  this  cannot  be  evaluated 
exactly  for  anything  other  than  the  simplest  curves,  it  lends  itself  naturally  to  asymptotic 
approximation.  Orthogonal  distance  regression  corresponds  to  a  first  order  approximation 
to  the  maximum  likelihood  estimator  in  this  model:  the  paper  also  derives  a  second  order 
approximation,  which  turns  out  to  be  a  simple  modification  of  the  least  squares  penalty 
that  includes  a  contribution  from  the  curvature  at  the  closest  point.  Analytical  and  nu¬ 
merical  examples  are  presented  to  demonstrate  the  improved  fits  possible  using  the  higher 
order  estimator,  along  with  an  extension  to  higher  dimensional  problems  involving  surface 
fitting. 
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2  Following  [10],  it  can  be  seen  that  points  on  the  curve  with  added  indepen¬ 

dently  and  identically  distributed  noise  are  more  likely  to  fall  in  the  concave 
region  containing  P'  than  in  the  convex  region  containing  P.  The  three  cir¬ 
cles  shown  with  centre  on  the  curve  represent  contours  of  constant  probability 
under  the  standard  stochastic.  Obviously,  more  points  on  the  curve  are  likely 
to  contribute  to  noisy  observations  at  P'  than  at  P.  Consequently,  an  obser¬ 
vation  is  more  likely  to  occur  at  P'  and  similarly  for  other  points  in  a  concave 
region  when  compared  with  those  in  a  convex  region .  9 

3  The  ellipse  with  a  —  1  and  b  =  0.2  is  sampled  uniformly  along  its  arc  length 
and  then  iid  zero-mean  Gaussian  noise  with  a  =  0.05  is  added  to  the  x  and 
y  coordinate  of  the  sampled  point  to  produce  the  5000  noisy  points  shown  in 
(a)  above.  The  contours  of  the  sum  of  squared  orthogonal  distances  (2)  and 
the  second-order  approximation  to  the  log-likelihood  function  (13)  are  shown 

in  (b)  and  (c)  respectively .  19 

4  The  ellipse  with  a  =  1  and  b  —  0.1  is  sampled  uniformly  along  its  arc  length 
and  then  iid  zero- mean  Gaussian  noise  with  a  =  0.05  is  added  to  the  x  and 
y  coordinate  of  the  sampled  point  to  produce  the  5000  noisy  points  shown  in 
(a)  above.  The  contours  of  the  sum  of  squared  orthogonal  distances  (2)  and 
the  second-order  approximation  to  the  log-likelihood  function  (13)  are  shown 

in  (b)  and  (c.)  respectively .  20 

A1  A  typical  ionogram.  Frequency  is  plotted  along  the  x  axis  and  signal  travel¬ 
time  along  the  y-axis.  The  grey-level  indicates  the  strength  of  the  received 
signal . 


25 


DSTO-RR-0242 


1  Introduction 

Many  problems  in  computer  vision  and  elsewhere  centre  on  finding  some  member  of 
a  family  of  curves  that  best  fits  a  given  set  of  points  in  a  plane  (or,  more  generally,  on 
finding  a  member  of  an  m-dimensional  family  of  surfaces  that  best  fits  a  set  of  points 
in  n-dimensional  space).  At  present  the  preferred  method  for  this  is  orthogonal  distance 
regression:  in  this  approach  the  best  fit  is  taken  to  be  the  curve  which  minimises  the 
sum  of  squares  of  distances  from  the  data  points  to  their  nearest  neighbours  on  the  curve. 
Orthogonal  distance  regression  is  well  known  to  be  computationally  difficult,  and  although 
some  good  algorithms  are  now  available  [8,  1],  it  can  still  be  very  hard  to  numerically 
determine  the  optimal  fit.  This  paper  is  not  concerned  with  efficient  computation,  however. 
Rather  it  explores  the  more  fundamental  problem  of  the  flaws  in  the  stochastic  model 
on  which  orthogonal  distance  regression  is  based,  and  develops  an  improved  model  and 
associated  distance  measures  that  produce  better  fits. 

The  paper  stemmed  from  the  authors’  practical  experience  in  fitting  curves  to  ionogram 
data  (see  [14,  15,  16,  17]  and  Appendix  A  for  a  brief  discussion  of  ionograms  and  their 
approximation),  while  the  approach  taken  was  inspired  by  some  perceptive  remarks  by 
Kanatani  on  the  distinction  between  geometric  and  statistical  distance  measures  in  fitting 
curves  to  data  [11,  12,  10,  p.  359].  Figure  1  shows  a  typical  trace  in  an  ionogram,  together 
with  an  idealized  fit  to  the  data.  It  was  found  that  curves  fitted  to  these  data  sets  by 
orthogonal  distance  regression  had  a  tendency  to  underfit  the  nose  of  the  trace  (i.e., 
not  to  extend  as  far  as  the  data).  Theoretical  confirmation  of  this  bias  in  the  fitting 
procedure  would  have  serious  implications  for  ionogram  modelling  due  to  the  importance 
of  accurately  locating  the  exact  extent  of  the  nose  in  applications  of  interest. 

The  tendency  of  orthogonal  distance  regression  to  underfit  corners  has  been  noted 
previously.  Kanatani’s  intuitive  explanation  for  this  is  that  the  sum  of  squared  distances 
is  primarily  a  geometric  measure  of  fit  and  does  not  account  for  all  the  stochastic  processes 
that  play  a  part  in  generating  real  data.  To  see  this,  consider  an  simple  example  in 
which  the  data  are  noisy  measurements  of  points  that  are  sampled  randomly  along  some 
underlying  curve  with  a  nose  like  that  in  figure  1.  Then  data  points  are  more  likely  to 
be  found  inside  the  nose  than  outside  it,  since  a  point  inside  the  nose  at  a  fixed  distance 
from  the  curve  will  be  closer  on  average  to  the  rest  of  the  curve  than  is  a  point  at  the 
same  distance  but  on  the  outside  of  the  curve.  This  bias  is  a  product  of  the  measurement 
process,  the  distribution  of  samples  on  the  true  curve,  and  the  curvature  of  the  nose:  the 
sum  of  squared  distances  does  not  take  all  of  these  issues  fully  into  account.  In  particular, 
a  curve  chosen  under  this  metric  is  likely  to  underfit  the  true  nose  since  there  are  likely  to 
be  more  data  points  inside  the  nose  than  outside. 

Kanatani’s  insight  prompted  us  to  develop  the  simple  but  reasonable  statistical  model 
presented  here  for  the  distribution  of  a  set  of  noisy  observations  as  a  function  of  the  un¬ 
derlying  curve.  We  then  derive  the  likelihood  function  for  observed  data  under  this  model: 
while  it  cannot  be  exactly  evaluated  for  anything  other  than  the  simplest  curves,  it  lends 
itself  naturally  to  asymptotic  approximations.  In  particular  we  show  that  under  this  model 
orthogonal  distance  regression  can  be  viewed  as  a  first  order  (i.e.,  linear)  approximation 
to  the  true  maximum  likelihood  estimator.  Moreover  a  second  order  approximation  gives 
an  estimator  that  is  a  simple  extension  of  orthogonal  regression  that  now  includes  con- 
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Figure  1:  An  ionogram  trace  overlayed  with  an  idealized  fit  using  the  curve  c  from  some 
family  C  of  curves. 


tributions  from  the  curvatures  at  the  closest  points  on  the  curve:  these  contributions 
compensate  for  the  tendency  already  noted  to  underfit  corners. 

Kanatani’s  argument  also  suggests  that  orthogonal  distance  regression  will  fail  to  accu¬ 
rately  reconstruct  globally  convex  curves,  although  in  this  case  the  argument  now  implies 
that  the  regression  will  overfit  the  true  curve.  To  see  this,  note  that  a  random  perturbation 
of  a  point  on  a  convex  curve  is  more  likely  to  fall  outside  the  curve  than  within  it.  There¬ 
fore,  under  the  simple  model  proposed  above,  more  data  points  are  likely  to  lie  outside  the 
curve  than  inside,  so  orthogonal  distance  regression  is  likely  to  produce  an  estimate  that 
lies  outside  the  actual  curve.  While  this  behaviour  does  not  seem  to  have  been  explicitly 
noted  in  the  literature  before,  it  is  certainly  apparent  in  the  examples  considered  here. 
Once  again,  however,  the  problem  is  addressed  by  a  proper  statistical  analysis.  In  the 
model  presented  here  the  likelihood  function  turns  out  to  include  a  teim  involving  the 
length  of  the  curve  being  fitted.  This  term  favours  shorter  curves,  counterbalancing  the 
asymmetry  in  the  distribution  of  observed  points. 

More  formally,  it  is  known  (see,  e.g.,  [6,  p.  247]),  if  not  widely  appreciated,  that 
orthogonal  distance  regression  may  be  both  biased  and  inconsistent.  By  this  we  mean 
that  the  average  of  fits  derived  from  repeated  finite  samples  of  a  given  underlying  cuive 
may  not  be  the  underlying  curve,  and  that  as  the  number  of  samples  tends  to  infinity 
while  the  noise  level  stays  the  same  the  fitted  curve  may  not  converge  to  the  true  curve. 
One  explanation  for  this  is  that  orthogonal  distance  regression  is  maximum  likelihood 
estimation  under  a  simple  statistical  model  in  which  the  location  of  the  samples  on  the 
true  curve  are  unknown  parameters.  Since  the  number  of  free  parameters  in  this  model 
grows  in  proportion  to  the  number  of  data  points,  none  of  the  classical  convergence  results 
for  maximum  likelihood  estimation  apply.  Another  way  of  viewing  the  model  presented 
here  is  that  it  extends  the  simple  model  in  a  way  that  allows  elimination  of  these  fiee 
parameters 

The  structure  of  the  paper  is  as  follows.  In  the  next  section  we  introduce  notation  and 
describe  the  standard  stochastic  model  of  data  formation:  the  observed  data  are  noisy 
perturbations  of  unknown  points  on  some  underlying  true  curve.  We  then  extend  this 
model  with  the  additional  assumption  that  the  unknown  points  are  uniformly  distributed 
along  the  true  curve.  Within  each  model  we  can  then  define  the  best  fitting  curve  as 
that  curve  which  maximises  the  likelihood  of  observing  the  given  data.  Under  Gaussian 
noise,  for  the  first  model  this  maximisation  is  just  orthogonal  distance  regression;  in  the 
second  model  the  likelihood  function  for  each  observation  takes  the  form  of  an  integral  of  a 
Gaussian  kernel  along  the  given  curve,  with  the  maximum  of  the  kernel  being  at  the  closest 
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point  on  the  curve  to  the  observation.  While  this  integral  is  intractable  in  general,  its  form 
immediately  suggests  approximation  by  an  asymptotic  expansion  about  the  closest  point. 

Section  3  derives  the  simplest  such  approximation  by  taking  a  locally  linear  approxima¬ 
tion  to  the  curve  at  the  closest  point  and  assuming  that:  the  data  points  are  independent 
samples;  the  variation  in  curve  length  is  insignificant;  and  the  radius  of  curvature  is  large 
compared  to  the  noise  level.  Under  this  approximation  the  maximisation  reduces  to  or¬ 
thogonal  distance  regression.  This  shows  maximum  likelihood  estimation  under  the  new 
model  is  a  credible  extension  of  orthogonal  distance  regression. 

Section  4  then  derives  an  improved  estimate  based  on  a  quadratic  approximation  of  the 
curve  at  the  closest  point,  and  shows  that  this  leads  to  the  addition  of  a  linear  correction 
to  each  squared  distance.  This  correction  term  is  product  of  the  distance  to  the  curve  and 
the  ratio  of  the  noise  variance  and  the  local  curvature.  Thus,  unlike  standard  orthogonal 
distance  regression,  the  penalty  function  requires  knowledge  of  the  variance:  this  is  usually 
unknown,  but  can  be  estimated  by  including  it  as  a  variable  in  the  likelihood  maximisation 
and  solving  for  it  explicitly. 

Section  5  elucidates  the  differences  between  the  fits  calculated  from  the  various  ap¬ 
proximations  derived  in  sections  3  &  4  for  the  special  case  of  fitting  a  circle  centred  on  the 
origin  to  observed  data.  In  this  case  analytic  expressions  can  be  derived  for  the  various 
maximum  likelihood  estimators,  along  with  expressions  for  the  bias  and  variance  in  these 
estimators.  The  bias  is  expressed  as  a  series  in  the  ratio  of  the  variance  to  the  square  of 
the  radius,  and  is  independent  of  the  number  of  observations.  Under  orthogonal  distance 
regression  the  bias  expansion  starts  with  a  first  order  term:  in  the  new  estimator  this  term 
is  eliminated  from  the  expansion.  We  also  determine  the  conditions  under  which  the  bias 
in  orthogonal  distance  regression  is  large  compared  to  the  variance,  and  so  under  which  it 
would  pay  to  adopt  the  more  accurate  estimator  proposed  here. 

Section  6  gives  some  more  realistic  illustrative  numerical  results  on  fitting  ellipses. 
These  show  the  tendency  of  orthogonal  distance  regression  to  overfit  convex  curves,  and 
that  the  second  order  estimator  presented  here  produces  a  more  accurate  fit.  Finally, 
after  conclusions  in  section  7,  some  supporting  material  is  presented  in  the  appendices. 
This  includes:  a  brief  description  of  the  autoscaling  problem  for  ionograms;  closed  form 
expressions  for  various  moment  integrals  used  in  the  main  body;  and  an  extension  of  the 
analysis  to  the  problem  of  fitting  an  d- 1  dimensional  surface  to  a  set  of  points  in  . 


2  The  Problem 

In  the  following  subsections,  we  firstly  examine  the  standard  model  and  show  how  it 
reduces  to  orthogonal  distance  regression,  and  secondly  we  extend  the  model  to  address 
the  failings  of  the  standard  model. 


2.1  The  Standard  Model 

Formally,  the  problem  of  interest  is  the  following: 
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Problem:  Given  a  set  of  noisy  observations  xi, . . .  ,  xjv  G  M2 ,  find  the  best  fitting  curve 

c  to  the  data  from  some  family  C  of  curves. 

As  noted  by  Kanatani,  a  necessary  precondition  for  any  well-defined  fitting  procedure  is 
a  model  for  data  formation  and  the  introduction  of  errors.  The  usual  model  is  [4,  12]: 

Standard  stochastic  model:  Each  x7l  is  an  independent  noisy  observation  of  a  point 
yn  on  the  curve  that  has  been  corrupted  by  additive  isotropic  Gaussian  noise  with 
variance  a 2. 


Note  that  this  model  has  two  distinct  components,  each  of  which  needs  to  be  deter¬ 
mined  in  the  fitting  process.  The  first  and  obvious  component  is  the  unknown  curve  c: 
this  will  usually  be  represented  as  some  function  c(0)  of  a  vector  6  =  ,0k)  of 

parameters.  The  second,  less  obvious,  component  is  the  nature  of  the  stochastic  processes 
determining  the  observed  distributions.  This  too  will  normally  be  restricted  to  some  known 
class  of  models  parameterised  by  a  vector  a  =  (cti,  . . .  ,  ol)-  The  fitting  process  must  then 
determine  the  parameters  6  and  <r. 

For  simplicity,  for  now  we  assume  the  noise  has  zero  mean  and  is  identically  and 
independently  distributed:  thus  it  can  be  characterised  solely  by  its  standard  deviation 
a.  In  many  practical  situations,  however,  the  noise  will  be  correlated;  for  example,  trace 
extraction  in  ionogram  autoscaling  and  edge  extraction  both  produce  connected  segments 
of  pixels  which  arc  then  to  be  fitted  by  lines  or  curves.  In  this  case  neighbouring  pixels 
on  the  segments  must  be  adjacent,  and  so  are  strongly  correlated:  how  this  correlation 
should  be  incorporated  into  the  model  will  be  considered  in  a  later  paper. 

Thus,  given  a  curve  c  and  noise  level  <7,  the  standard  model  defines  the  conditional 
probability  p(x|y;c,  <x)  that  a  point  x  is  a  noisy  observation  of  a  point  y  on  c  as 


1 


-llx-yll2/2(T2 


subject  to  the  constraint  y  G  c.  This  gives  the  following  joint  distribution  for  data  and 
model: 


N 

p(x i,...  ,Xjv|yi,...  ,yjv;c,cr)  =  JJp(xn|yn;c,cr). 

i— 1 


Under  this  model  the  curve  c,  noise  level  a  and  unknown  true  points  yi,...  ,  y n  aie 
all  parameters  describing  a  distribution  that  are  to  be  inferred  from  the  available  data. 
Maximum  likelihood  estimation  has  them  being  chosen  to  minimise  the  log-likelihood 
function 

N 

l (xx , . . .  ,xjv  :  yi, . . .  ,y n,c,o)  =  -  lnp(x„|y„;  c,  a)  (1) 

11  — \ 


giving  the  problem 


min 
yi . y 


1 

2^ 


N 


iix»  ~ y?i 


+  2N  lncr 


subject  to  yi,...  ,yv  Gc. 


(2) 
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It  is  clear  that  for  fixed  c  the  optimal  choice  for  y„  is  yn  =  <,  where  x*t  is  the  closest 
point  on  c  to  xn.  Moreover,  the  yn  are  independent  of  a,  and,  if  a  is  unknown,  optnmsmg 
(2)  with  respect  to  a  gives  that  its  maximum  likelihood  estimate  is 

0*  =  — -Yd2n.  (3) 

2  N^n 

n=l 

Here  dn  =  ||xn  -  x*  ||  is  the  distance  from  xn  to  c.  Substituting  this  expression  for  a  back 
into  (2)  and  exponentiating  reduces  it  to  standard  orthogonal  distance  regression. 

While  this  analysis  provides  a  statistical  foundation  for  orthogonal  distance  regression 
it  also  highlights  the  method’s  drawbacks.  In  particular,  it  explicitly  requires  estimation  o 
the  unknown  true  sample  points  yn  as  well  as  c  and  a.  In  practice  the  points  yn  are  of  httle 
interest  in  themselves,  and  a  model  for  the  reduced  conditional  density  P(xx, .  .  ,  x*  c,  a) 
would  be  preferable  to  one  for  p(xx, . .  •  ,xjv|yi, .  •  •  ,  Yn;  c,  a).  The  y„  can  only  be  e  im 
nated  within  the  standard  model,  however,  by  strong  assumptions  such  as  the  noise  level 
being  small  compared  to  the  curvature  (see,  e.g.,  [13,  12]).  We  therefore  propose  an 

alternative  model. 


2.2  The  Extended  Model 


Ideally  curve  estimation  would  be  based  on  a  complete  statistical  analysis  that  started 
from  a  model  for  the  full  joint  distribution  p(x,y;c,a).  This  can  be  factored  as 

p(x,  y;  c,  or)  =  p(x|y;  c,  a)p(y;  c,  a). 

The  standard  model  stops  here,  providing  an  expression  for  p(x|y;c,<r)  but  making  no 
assumptions  about  the  form  of  p(y;  c,  a).  We  now  extend  the  model  by  noting  that,  under 
the  reasonable  assumption  that  the  curve  generation,  curve  sampling  and  noise  piocesses 
are  independent,  p(y;  c,  a)  can  be  further  factored  as 

p(y;c,cr)  =p(y|c)p(c)p(cr). 

We  do  not  propose  a  complete  model  here,  as  the  prior  distributions  p(c)  and  p(o)  are 
very  problem  dependent  and  usually  impossible  to  specify  m  advance.  We  can,  howevei, 
propose  a  reasonable  universal  model  for  the  sampling  process  p(y|c)  that  generates  the 
points  observed  on  the  curve,  in  particular  that  it  is  uniform  in  arclength.  Thus 


p(y|c)  = 


ds  ds 


where  ds  is  the  element  of  arc  length  and  |c|  is  the  length  of  c. 

With  this  model  the  unknown  true  sample  points  can  now  be  removed  from  the  analysis 
by  noting  that 


p(x|c,cr)  =  j  p(x|y,c,<r)p(y|c)dy 

.  _  [eHI*-y(s)ll2/2a2ds 

a2\c\  Jc 
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where  the  points  y(.s)  along  the  curve  have  been  parameterized  by  arc  length  s.  This  gives 
an  expression  for  the  conditional  probability 

N 

p(x  1,...  ,xjv|c,cr)  =  JJp(x„|c,cr)  (5) 

i— 1 

as  a  product  of  integrals  that  no  longer  involve  the  y„. 

We  are  now  in  a  position  to  define  the  maximum  likelihood  estimate  of  the  best  fitting 
curve  c  to  the  data.  As  usual  we  define  the  log-likelihood  function 

f(x„  :  c,  a)  =  -  \np(xn\c,a). 

Therefore  the  maximum  likelihood  estimate  of  the  true  underlying  curve  c  is  the  curve  c 
that  minimizes 

N 

l{x  1,...  ,xAr  :  c,cr)  =  -  ^  lnp(x„  |c,  ct).  (6) 

71  =  1 

It  now  remains  to  evaluate  the  terms  in  (6).  In  practice  it  will  be  impossible  to  get 
closed  form  expressions  for  l(xn  :  c,cr)  for  all  but  the  simplest  of  curves,  so  we  shall  have 
to  resort  to  asymptotic  expansions  of  the  integral  in  (4).  But  before  starting  this,  we  first 
wish  to  make  a  number  of  further  comments  on  the  model  as  proposed  so  far. 

First,  note  that  the  above  analysis  shows  that  the  a  posteriori  probability  of  a  point  y 
on  a  curve  c  generating  a  given  observation  x  is 

,  ,  X  p(x|y,a)p(y|c) 

p(y  x,c,  cr)  =  - : — ; - 7 - 

’  ’  '  p(x|c, a) 

e-||x-y||2/2<x2 

=  Jc  e“llx_y(s)H2/2ff2ds ' 

Second,  in  many  cases  the  admissable  curves  c  will  have  infinite  length,  for  example 
when  C  is  a  family  of  straight  lines.  I11  this  case  we  shall  model  p(y  |c)  by  an  non-inf ormative 
prior  [3].  I11  this  model  we  assume  all  points  on  c  are  equally  likely  and  all  curves  in  C 
have  the  same  length,  and  we  interpret  these  statements  to  mean  that  the  factor  |c|  can 
simply  be  removed  from  (6). 

Third,  in  principle  the  ideal  quantity  to  compute  is  the  a  posteriori  density 

p(c,a ) 

p(c,a|Xl,...  ,XJV)  =P(X1,...  ,xN\c,a)p(xi  "  x-~y 

Knowledge  of  p(c,a|xi, . . .  ,xN)  would  give  a  much  better  overall  understanding  of  the 
problem,  and  would  allow  construction  of  Bayesian  estimators  that  minimised  appropriate 
loss  functions.  Note  that  in  practice  the  curve  generation  and  noise  generation  processes 
will  most  likely  be  independent,  so  that  p(c,a )  =  p(c)p(cr). 

Unfortunately  computation  of  p(c ,  tr|xi, . . .  ,xa 0  requires  knowledge  of  the  prior  dis¬ 
tributions  p(xi, . . .  ,xjv),p(c)  and  p(<r),  as  well  as  p(xi, . . .  ,xyv| c,a).  In  some  cases  this 
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may  be  available.  For  example,  often  there  will  be  no  preferred  ongm  m  the  problem  o 
that  the  family  C  inclndes  all  translates  of  itself.  In  this  case  p(x,, . . .  ,x„)  will  be  uniform 
regardless  ofp(c)  and  p(a).  On  the  other  hand,  if,  as  is  often  the  case  the  observations 
are  first  scaled  and  translated  to  the  unit  square  before  fitting,  then  p(xi,  •••.*«)  wil 
not  be  uniform  but  will  reflect  the  prior  distribution  of  curves  p(c). 

More  importantly,  it  is  not  at  all  clear  what  form  a  reasonable  prior  distribution  p(c) 
on  the  space  of  all  curves  might  take,  or  how  to  define  it.  Possibilities  exist  such  as  Weiner 
measures  on  spaces  of  curves,  but  a  discussion  of  them  is  well  beyond  the  scope  of  tl 

paper. 

Finally,  as  noted  before,  in  practice  in  most  fitting  problems  the  class  of  curves  C 
consists  of  a  family  of  shapes  c{0)  parameterized  by  the  vector  6  -  (6U...  M-  n  n 
case  maximum  likelihood  estimation  now  becomes  maximisation  of  p(xi,. . .  ,xN  W)  ovei 
0,  while  the  expression  for  the  posterior  distribution  is 

p(0) 

p(0|Xl,...  ,xjv)=p(x  1,...  ,xN|0)^-— 

Knowledge  of  the  posterior  now  requires  knowledge  of  the  prior  distribution  on  the  pa¬ 
rameter  space.  This  can  be  even  more  problematical  than  knowledge  of  the  prior  on  some 
space  of  curves:  with  a  poor  (or  even  a  natural)  parameterization,  the  natural  measure  on 
the  parameter  space  may  not  be  very  closely  related  to  the  natural  measure  on  the  space 
of  actual  curves.  Nevertheless  the  posterior  distribution  p(0|x)  would  still  be  the  prefeired 
distribution  to  work  with  if  any  further  analysis  is  to  be  carried  out,  such  as  hypothesis 
testing  to  decide  if  curves  can  be  merged. 


3  First  order  approximation 


To  compute  the  likelihood  function  we  need  to  evaluate  integrals  of  the  form 

1  ^  0-Hxn-y(s)H2/2lj2(iS. 


p(x„|c,cr)  = 


(7) 


2ixo2\c 

In  practice  this  will  be  impossible  to  calculate  exactly  for  all  but  the  simplest  of  curves. 
The  Gaussian  kernel,  however,  naturally  suggests  some  form  of  asymptotic  approximation 
along  the  general  lines  of  Laplace’s  method  [5]  around  the  maximum  of  the  mtegran  . 

This  occurs  at  the  point 

x*  =  argmin{||x„  -  y(s)(|2  :  y(s)  €  c} 

i.e.,  x*  is  the  closest  point  to  xn  on  c.  We  shall  denote  the  distance  ||xn  -  x*  ||  by  dn. 

We  begin  by  constructing  a  simple  linear  approximation  to  the  integral  in  (7).  we  shall 
show  that  under  this  maximum  likelihood  estimation  again  reduces  to  orthogonal  distance 
regression.  Let  n„  be  the  normal  to  c  at  <,  and  let  Tn  _  {y  :  nn(y  x  )  J  ® 
the  tangent  line  to  c.  Thus,  if  the  curve  is  smooth  and  varies  very  little  over  the  cen Lai 
support  of  the  Gaussian  kernel  {i.e.,  the  radius  of  curvature  is  large  compared  to  a),  then 
we  may  approximate  the  integral  in  (7)  by 

(x  ,c  a)  „  J_J_  [  e-^-^2/2a2ds. 

PlXn'  ’  j  27TCT2  |C|  Jy£Tn 
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Since  x*  is  the  closest,  point  to  x„  on  c,  x„  —  x*  must  be  perpendicular  to  Tn,  i.e.,  xn  —  x* 
must  be  proportional  to  n„.  Therefore  by  Pythagoras’  theorem  for  y  €  Tn, 

l|x„  -  y||2  =  ||x„  -  x*  ||2  +  ||x*  -  y||2, 


from  which  it  follows  that 

p(xn\c,a) 


1  e-d?,/2  a2 

2na2  |c| 

1  e-dl/2a2 


/, 


-IK-yll2/2or2^s 


\/2na  |c| 

Thus  the  log-likelihood  function  reduces  to 


"  dl 


l{xu. . .  ,xyv  :  c,a)  ~  E  +  N  |  In  |c|  +  In  V2na 


n— 1 


(8) 


This  loss  function  is  identical  with  that  used  in  orthogonal  distance  regression,  except 
for  the  additional  contributions  due  to  the  variable  curve  length  and  the  unknown  noise 
level.  If  both  of  these  were  known  and  fixed,  then  the  linear  approximation  to  maximum 
likelihood  estimation  would  indeed  reduce  to  orthogonal  distance  regression.  In  practice, 
however,  both  terms  will  vary,  and,  as  shall  be  seen  in  section  5,  their  inclusion  can  make 
a  significant  difference  to  the  results  produced  by  straight  orthogonal  distance  regression. 

In  further  analysing  the  contributions  from  the  terms  involving  the  variance  and  curve 
length  to  the  maximum  likelihood  estimator,  note  first  that  the  variance  must  be  known 
or  estimated  if  the  curve  length  contribution  is  included.  This  follows  since  the  relative 
weight  of  the  contributions  from  the  sum  of  squared  differences  and  log  of  curve  length  is 
moderated  by  the  variance.  Fortunately,  if  the  variance  is  unknown  it  can  be  eliminated 
by  direct  substitution.  To  see  this,  note  that  at  a  minimum  of  l,  dl/da  =  0  must  hold 
regardless  of  the  form  of  c.  This  implies 


a 


2 


1 

N 


N 


71=1 


Substituting  this  back  into  (8),  disregarding  constants,  exponentiating  and  then  squaring 
shows  that  the  first  order  approximate  maximum  likelihood  estimate  of  c  is  the  minimiser 
of 

N 

ei(xi, . . .  ,x/y  :  c)  =  |c|2  ^d2.  (9) 

71  =  1 

While  this  minimisation  is  still  very  close  to  orthogonal  distance  regression,  curve  length 
now  obviously  plays  a  significant  role. 

The  effect  of  the  curve  length  contribution  has  a  reasonable  interpretation:  all  other 
things  being  equal  (i.e.,  if  the  closest  point  remains  unchanged  and  the  first  order  asymp¬ 
totic  approximation  is  still  valid),  a  shorter  curve  has  a  smaller  range  of  y  values  that 
could  have  given  rise  to  the  observation  and  therefore  an  increased  probability  that  a 
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Koure  l-  Following  [10],  it  can  be  seen  that  points  on  the  curve  with  added  independently 
and  identically  distributed  noise  are  more  likely  to  fall  in  the  concave  region  containing 
P'  than  in  the  convex  region  containing  P.  The  three  circles  shown  with  centre  on  the 
curve  represent  contours  of  constant  probability  under  the  standard  stochastic.  Obviously, 
mom  points  on  the  curve  are  likely  to  contribute  to  noisy  observations  at  P  than  at  P 
Consequently,  an  observation  is  more  likely  to  occur  at  P'  and  similarly  for  other  points 
in  a  concave  region  when  compared  with  those  in  a  convex  region. 


neighbourhood  of  the  closest  point  was  indeed  the  source  of  the  observation  And  while 
ei  XJV  :  c)  can  be  decreased  to  zero  if  |c|  |  0,  this  behaviour  violates  the  assump¬ 

tions  under  which  (8)  and  (9)  were  derived,  as  the  curvature  will  no  longer  be  large  with 
respect  to  the  noise  level  at  any  point  xn. 

Strictly  speaking,  a  full  one-dimensional  approximation  would  see  the  terms  involving 
|C|  disappear  from  (8)  and  (9):  since  tangent  lines  have  infinite  length  and  the  length  of 
Lll  tangents  is  the  same,  as  noted  in  the  previous  section  it  would  be  reasonable  to  use  a 
non-informative  prior.  This  has  its  attractions,  especially  when  |c|  is  difficult  to  compute 
and  the  curves  under  consideration  have  infinite  length  and  no  well  defined  boundary. 
It  conflicts  with  strict  maximum  likelihood  esimation,  however,  and,  as  section  5  shows, 
terms  involving  |c|  can  make  significant  contributions. 

Finally  we  note  again  that  in  practice  i(*r . XN  :  c,<r)  will  usually  be  represented 

in  parameters  form  as  !(*„.. .  ,x«  :  0,o),  and  the  problem  actually  faced  w.ll  be  to 
minimise  (8)  or  (9)  as  functions  of  0. 


4  Second  Order  Approximation 

The  problem  with  the  first  order  model  in  that  it  does  not  account  for  the  effect  that 
local  curvature  is  likely  to  have  on  the  distribution  of  observed  points.  For  example, 
consider  the  following  situation  of  two  points  P  and  P'  that  are  equidistant  to  an  on 
either  side  of  the  nose  of  a  curve,  and  such  that  the  nose  is  the  closest  point  on  the  curve 
to  both  points  (figure  2,  [10]).  The  positive  curvature  means  that  point  P  is  more  likely 
to  be  generated  by  the  noise  process  than  point  P  since  P  is  closer  to  the  legs  of  the 
curve.  Thus  while  each  point  will  make  an  equal  but  opposite  contribution  to  locating  t  le 
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both  points  will  have  an  equal  but  opposite  effect  in  positioning  the  curve,  u  ^  q 

and  the  curvature  of  the  nose:  it  is  not  properly  accounted  for  in  minimising 
squared  distances. 

A  similar  argument  suggests  that  orthogonal  distance  regression  *0  gij  a  b,a*id 
it,.  Moreover  points  falling  outside  aie  more  like  y  to  e  i  f  c  tj  c  observed 

that  this  is  indeed  the  case. 

The  above  argutnonts,  that  orthogonal  distance  regression  based  o„ the  linear  approx- 
i.nation  may  produce  biased  fits,  are  supported  by  our  observat »  tl.at^  ^ 

using  orthogonal  distance  regression  can  be  biased  mound  the  i  _  log-likelihood 

ionograms  [14,  15,  10,  17].  Therefore  we  seek  a  better  approximation  to  the  log 

function  that  will  take  into  account  the  effect  of  local  curvatuie. 


4.1 


Deriving  the  second  order  approximation 


Let  c*  =  <(r 


be  the  best  approximating  circle  to  c  at  x*  ,  and  suppose  that  c. 

Again  x„  -  xn  is 


has  centre  z, 


and  radius  rn  (rn  is  the  inverse  of  the  curvature  of  c  at  x 


U1U4  ill  V  oroc  wx  - - 

translate  and  rotate  the  plane  so  that  the  global  coordinate  system  cornu  cs  wnao 
coordinate  system  in  which 

Z  n  =  (0,0), 
xn  =  ( rn  —  dn,  0), 

X*  =  {l'n ,  0). 

We  shall  adopt  the  convention  that  d„  is  signed.  This  will  indicate  wlmther ^  cohere 
or  convex  with  repect  to  x„:  if  *.  is  posit, ve  the  x„  res  “ms.de  ^  ^ 

[fir'ftl^signS'^t  n'^'Z *  indicate'tlie  local  geometry,  and  indeed 

integrals  that  appear  below.) 

With  this  approximation  we  have 


p(x„  |c,  o) 


L-Lf 

TCT2  |c|  Jc’n 
2na2  |c|  f0 


2ira 

1 


-y(*)ll 2/2a  ds 


-(Zn+rnicasO’™'6))^  7'd0 
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and 

||xn  -  (zn  +  rn(cos6>,sin0))||2  =  \\{rn  ~  dn, 0)  -  (r„  cos6>,  rn  sm0)|| 

=  {rn-dn-rncosd)2  +  rlsin29 
=  2rn{rn-dn){l-cos0)  +  dl. 

The  integral  can  be  be  evaluated  exactly  as  a  Bessel  function,  giving 

p(xn\c,a)  ~  2^2 
_  l_Tn 

”  a2  |c| 


rn  -dl 

- — rC 


£/2 (T2  e-rn(rn-fin)/(r2 

JO 


dd 


-dH2a*e-rn{rn-dn)la*lo  jn)^  . 


(10) 


We  now  look  at  an  asymptotic  expansion  of  p(xn|c,a)  under  the  assumptions 

«  i  dr 

>  1, 


a 


£»1, 

dn. 


and 


a 


—  H  ~  1. 


These  assumptions  are  consistent  with  the  underlying  model:  we  expect  that  the  distance 
beZel  a  typical  observation  x„  and  the  closest  point  on  the  curve  to  ,t  w.ll  be  of  order 
a.  Thus,  after  defining  p  =  rn/a ,  we  seek  an  asymptotic  expansion  o 

^(p)=pe"^-M)/o(p(p-A*)) 


for  large  p.  To  order  p  2  we  have 

M  ~ 


pp(p-v) 


\/27r p(p  -  p) 


1  1 
1  +  x 


8p(p~p)J 


giving 


—  In  </>(p)  ~  ln  ^ 

~  In  V27T  -  - 


fi-el 

-ln 

L  P- 

1  + 


1 


1 


1 

1 

2 

P 

2 

w 

8p2  (1  —  M/P)  J 

1 

8?  ’ 


(11) 


and  therefore  that 

—  lnp(xn|c,  a) 


dl 


ln\/27r cr  +  In  |c|  +  ^ 


dn 

2rn 


8r2 


+  2  dr 


(12) 


It  is  not  obvious  that  the  coefficients  of  the  second  order  terms  a2/r2  and  d2/r2 
A  .1  ,1 ko  uHoroH  Kv  a  more  ac.CUI 


m 


It  is  not  obvious  tnai  me  wcmoioirto  ^  ox -  ■ 

the  above  equation  are  correct  as  is,  or  whether  they  would  be  altered  by  a  more  accurate 
approximation  to  c  that  involved  higher  derivatives.  Therefore  we  shall  only  mclude  the 

first  order  term  to  give 


/(x1? .  •  •  ,XJ V  :  c,cr) 


N 

£ 

n= 1 


2ct2 


dn 
2  rn 


+  ./V  In  |c|  +  ln  \/27T(tJ 


(13) 
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(However  sec  the  end  of  the  next  subsection  for  a  discussion  of  an  alternative  representation 
of  the  log-likelihood  function  that  includes  some  second  order  terms). 

Equation  (13)  has  one  significant  difference  from  (8):  the  linear  correction  term  dn/2rn. 
This  will  act  to  counter  the  problem  of  underfitting  noses  noted  earlier.  For  if  xn  lies  inside 
the  curve  then  dn  is  positive,  so  the  decrease  in  the  term  dn/rn  will  offset  to  some  degree 
the  increase  in  d2J2a2  as  the  curve  is  moved  away  from  x„ .  In  contrast,  if  x„  lies  outside 
the  curve  then  both  d2/2a2  and  dn/rn  will  be  positive,  and  so  will  attract  the  curve 
towards  x„. 


4.2  Alternative  forms  for  the  second  order  approximation 

We  now  present  some  alternative  formulations  of  (13)  that  may  be  more  appropriate 
for  computational  purposes,  or  that  give  further  insight  into  the  underlying  concepts. 
First,  an  unknown  variance  can  again  be  eliminated  from  (13)  by  direct  substitution 
of  the  maximum  likelihood  estimate  of  c.  Since  a  does  not  appear  in  the  linear  term, 
the  condition  dl/da  =  0  implies  the  maximum  likelihood  estimate  of  the  variance  is 
a2  =  jj  ZnLi  dn-  Substituting  this  back  into  (13),  disregarding  constants,  exponentiating 
and  then  squaring  gives  that  the  second  order  approximate  maximum  likelihood  estimate 
of  c  is  the  mini  miser  of 


C2(xi ,  .  .  .  ,XA r  :  c) 


•  exp 


1  N  A 
1  yv  Un 

N^~n 

71  =  1 


(14) 


Next,  it  is  worth  noting  that  slightly  different  versions  of  (13)  and  (14)  can  be  derived 
by  going  back  to  (11)  and  retaining  the  term  ln(l  —  p/p)  as  is,  rather  than  replacing  it  by 
a  first  order  approximation.  This  gives 


N 

l(x  1,...  ,XN  :  c,  a )  ~  ^ 

71  =  1 


d2  1 

#+2ln 


+  N  [in 


+  In  \/2tw  , 


and 


e2(xi,...  ,xAr  :  c)  =  |c|2-  (  ^  XX, 

\  71  =  1 

Nevertheless,  while  these  equations  may  be  more  accurate  approximations  to  the  true  log- 
likelihood  and  associated  functions  in  principle,  in  practice  (13)  and  (14)  arc  equivalent  to 
first  order  and  are  also  easier  to  minimise  numerically  since  they  allow  dn/rn  to  be  greater 
than  1  without  giving  rise  to  singularities  or  similar  difficulties. 

Finally  it  is  instructive  to  explore  the  form  of  (13)  a  little  further.  In  particular,  it  is 
tempting  to  complete  the  square  and  rewrite  the  equation  as 


N 

n 

Ln=l 


N 

l{xi,...  ,xN  :  c,  a)  ~ 

71=1 


(dn  -  <72/2rn)2 
2cr2 


+  N 


In  |c|  +  In  \/2iuj 


(15) 
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This  can  then  be  interpreted  as  the  log-likelihood  function  of  a  normal  distribution  approx¬ 
imation  to  p(xn\c,cr).  This  interpretation,  however,  would  be  wrong  and  misleading,  as 
well  as  not  agreeing  to  second  order  with  (12).  To  understand  why  it  suffices  to  note  that, 
to  first  order,  E[dn]  ~  —tj2/2rn:  this  follows  from  setting  pn  —  rn  —  dn  and  appealing  to 
(B2)  of  Appendix  B.  (Here  E[x\  denotes,  as  usual  the  expectation  of  the  random  variable 
x.)  Thus  the  term  a1  /2rn  in  (15)  cannot  be  interpreted  as  the  mean  of  p(x7l|c,  a)  as  it 
has  the  wrong  sign.  Rather  p(xn|c,  a)  is  a  skewed  distribution  and  cr2 /2rn  is  a  first  order 
approximation  to  its  maximum  (or  mode):  as  the  mode  is  the  relevant  statistic  in  maxi¬ 
mum  likelihood  estimation,  the  approximation  has  instead  constructed  a  local  quadratic 
approximation  to  the  likelihood  function  about  the  mode  rather  than  the  mean.  This  con¬ 
firms  and  clarifies  both  the  observations  made  in  the  introduction:  that,  while  on  average 
an  observation  is  more  likely  to  be  made  outside  the  curve  than  inside,  the  most  likely 
single  point  for  any  observation  is  just  inside  the  curve. 


5  Comparison  of  approximate  maximum 
likelihood  estimates  of  circles 

To  get  a  feel  for  the  differences  between  the  solutions  of  each  of  the  various  approximate 
maximum  likelihood  estimators  introduced  above,  we  apply  them  here  to  the  particular 
problem  of  fitting  a  circle  to  noisy  data.  Let  Co  be  a  circle  of  radius  ro  centred  on  the 
origin,  and  let  {xn}^=1  be  a  collection  of  noisy  observations  of  Co  generated  according  to 
the  basic  model  under  noise  level  a.  We  seek  a  circle  c  of  radius  r  also  centred  on  the 
origin  that  best  fits  this  data:  i.  e. ,  we  seek  an  estimate  r*  of  r0. 

To  achieve  this,  we  first  represent  the  data  in  polar  coordinates  as  xn  =  (pn  cos  9n,  pn  sin  0n). 
Then  substituting  ro  —  pn  for  dn  in  (10)  shows  that  the  distribution  of  the  xn  can  be  sum¬ 
marised  in  the  density  function 

p(xn|c0,  a)  =  — 1-^ e-(ro-^)2/2^2e-roPn/<T2/0  pn  dpn  ddn.  (16) 

We  now  use  this  distribution  to  calculate  the  estimates  that  will  be  returned  by  each  of 
the  approximations  considered  so  far. 


5.1  Orthogonal  distance  regression 

The  squared  distance  of  xn  to  any  circle  c  of  radius  r  centred  on  the  origin  is  simply 
( r  —  pn)2 ■  Thus  the  estimate  of  ro  returned  by  orthogonal  distance  regression  will  be  the 
solution  r*  of 

N 

mrin  -  ^)2- 
*  71  —  1 

Thus 

=  <17> 

71=1 
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Since  the  pn  are  random  variables,  so  also  is  r*.  In  particular, 

i°ilS)pdp 

p~e  ''  '  ""  * Iq  dp. 


E[r*]  =  E[pn]  =  ^  r  pe-^^e-^h 

°  Jo 

p-rH'la2  ro o 


Oc 


/ 


2  -p2/2«72 


This  integral  is  evaluated  in  Appendix  B:  under  the  assumption  that  ro/cr  >  1  (B2)  gives 


E[r*}  ~  r0 


i+K^)  +Kt) 


(18) 


Equation  (18)  shows  that  straight  orthogonal  distance  regression  will  overestimate  the 
true  radius.  This  is  not  unexpected:  as  noted  previously,  since  the  circle  is  a  convex  curve 
a  random  perturbation  of  any  point  on  it  is  more  likely  to  fall  outside  it  than  inside. 

In  a  similar  vein  the  variance  of  the  orthogonal  distance  regression  estimate  is 

r  *i  var [pn]  E[pl)  -  E[pn)2 
*»■[<■  i  =  -jj-  =  — jj - ■ 

Equation  (B3)  of  Appendix  B  gives  E[p2t]  —  r2  +  2cr2.  Therefore  retaining  all  terms  to 
order  (o/r q)4  in  (B2)  gives 


var[r*]  ~ 


_1_ 

N 


7-q  +  2a2 


r 2 
r0 


N 


(19) 


The  variance  estimate  shows  that  if  the  bias  in  the  orthogonal  distance  regression 
estimate  for  ro  predicted  by  (18)  is  to  be  less  than  the  estimate’s  variance  (and  thus  not 
be  the  dominant  term  in  the  error  in  the  estimate),  then  to  first  order  N  must  satisfy 

Op a  a 2 

y/N  "  2  r0 


so  that 


2apr0 
a 

Here  av  is  the  confidence  level  that  ensures  that  the  estimate  r*  lies  in  the  interval  [ro  + 
o2  j (2/-o)  -  (Xpo/sfN ,  7-o  +  ct2/(2? -0)  -  avo/\fN]  with  probability  p. 
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5.2  First  order  estimate  with  known  variance 


We  now  present  the  estimates  r*  that  will  be  returned  by  the  various  approximations  to 
the  maximum  likelihood  estimator  derived  in  previous  sections.  Given  a  finite  number  N 
of  observations,  all  of  these  estimates  will  be  random  variables.  They  are,  however,  now 
solutions  of  quadratics  or  higher  order  equations,  and  so  have  more  complicated  forms 
than  (17).  Therefore  we  shall  not  attempt  to  determine  their  exact  means  and  variances 
here.  Instead  we  shall  simply  determine  the  asymptotic  values  of  the  estimators  as  N  f  oo 
and  note  that  the  associated  variances  will  again  asymptote  to  a1  jN .  This  analysis  will 
be  sufficient  to  illustrate  the  differences  between  the  various  estimators. 


We  start  by  minimising  (8)  under  the  assumption  the  noise  level  a  is  known.  Dividing 
by  N,  taking  the  limit  as  N  t  oo  in  (8)  and  disregarding  constants  gives  a  corresponding 
estimate  of  ro  as  the  solution  r*  of 


min 

r 


r2  -  2 rE[pn)  +  E[p2n] 
2  a2 


+  In  r. 


It  is  straightforward  to  show  that 

...  %»]  +  (E\pn\2  -  4<t;)1/2 

2 

Making  the  usual  assumption  that  tq/o  3>  1  and  substituting  in  the  asymptotic  expansion 
for  E\pn\  given  in  (18)  shows  that  to  fourth  order 


r*  ~  ro  1 


1 

2 


3  fa 


4' 


8 


(21) 


Thus,  if  the  variance  is  known  and  the  proposed  stochastic  process  for  generating  ob¬ 
servations  is  correct,  including  a  contribution  from  curve  length  in  orthogonal  distance 
regression  will,  to  second  order,  see  the  resulting  estimate  undershoot  the  true  radius  by 
as  much  as  the  original  orthogonal  distance  regression  estimate  overshot  the  mark. 


5.3  First  order  estimate  with  unknown  variance 


If  the  noise  variance  is  unknown  then  (9)  can  be  used.  For  circle  fitting,  taking  the 
limit  as  N  t  oo  in  (9)  and  disregarding  constants  gives 

min  r2(r2  -  2 rE[pn]  +  E[p2n}). 
r 

Again  it  is  straightforward  to  show  that 


r  — 


3 E[pn]  +  (9 E[Pn}2  -  8 E[p2n))V2 


Under  the  assumption  that  a/ro  1,  using  (B3)  and  (18)  gives  that  to  fourth  order 


r  ~  ro 


2  \r0J  8  \r0 


(22) 


Thus,  if  the  variance  is  unknown,  substituting  the  maximum  likelihood  variance  estimate 
gives  an  estimate  that  agrees  to  first  order  with  the  estimate  that  uses  the  known  variance. 
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5.4  Second  order  estimates  with  known  variance 

Wo  now  derive  the  estimate  of  r*  given  by  a  second  order  approximation  when  the 
variance  is  known.  Dividing  by  N  in  (13),  then  taking  the  limit  as  IV  |  oo  and  disregarding 
constants  gives  a  corresponding  estimate  of  r$  as  the  solution  r*  of 

min  r*  ~  2'  Elpj  t.Ml  -  !— g M  + 
r  2 az  2r 

Tims  r*  must  be  a  root  of 

^52=1=0.  (23) 

We  next  make  the  same  assumptions  and  substitutions  for  E[pn]  as  in  the  corresponding 
first  order  case.  This  is  followed  by  expanding  r  as  a  series  in  powers  of  (a/ro)2  with 
unknown  coefficients,  substituting  this  series  into  (23)  and  then  choosing  the  unknown 
coefficients  to  ensure  that  the  first  few  terms  in  the  result  are  identically  zero.  This  gives 
that  to  fourth  order 


r*  ~  to 


(24) 


Thus,  as  expected,  using  the  second  order  approximation  to  the  log-likelihood  function 
eliminates  the  leading  order  term  in  the  bias  of  the  estimate  of  the  true  radius  7'o- 


5.5  Second  order  estimate  with  unknown  variance 

If  the  noise  variance  is  unknown  then  (14)  can  be  used.  For  circle  fitting,  taking  the 
limit  as  ./V  f  oo  in  (14)  and  disregarding  constants  gives  the  estimate  of  r o  is  the  solution 
r*  of 


min  r2(r2  -  2rE[pn]  +  E[pl\)eE[pn]/r . 

r 

Therefore  r*  must  be  a  root  of 

4 r3  -  6 r2E[pn]  +  2 E[p2x\  -  E[pn}(r2  -  2 rE[pn]  +  E[p2,])  =  0. 


We  again  make  the  same  assumptions  and  substitutions  as  in  the  first  order  estimate  with 
unknown  variance  and  the  second  order  estimate  above.  This  gives  that  to  fourth  order 


r*  ~  r0 


(25) 
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5.6  Conclusions  on  circle  fitting 

We  pause  here  to  summarise  the  results  presented  in  this  section.  The  analysis  clearly 
shows  that,  under  a  fixed  noise  level,  orthogonal  distance  regression  gives  a  biased  estimate 
of  the  true  radius.  Moreover  this  estimate  is  inconsistent,  in  that  the  bias  does  not  go 
to  zero  even  if  the  number  of  data  points  N  f  oo.  Moreover  the  nature  of  the  bias  is 
such  that  orthogonal  distance  regression  will  tend  to  overfit  the  curve:  i.e.,  it  will  return 
an  estimate  r*  of  the  true  radius  ro  with  r*  >  r o-  This  bias  is  likely  to  be  dominated, 
however,  by  the  variance  inherent  in  the  fitting  of  noisy  data  unless  N  exceeds  the  inverse 
of  the  relative  bias. 

Including  curve  length  in  orthogonal  distance  regression  over-compensates  for  the  ten¬ 
dency  to  overfit.  The  leading  order  term  in  the  bias  expansion  now  has  the  same  magnitude 
but  opposite  sign  as  the  corresponding  term  in  the  bias  expansion  for  the  straight  orthog¬ 
onal  distance  regression  estimate.  Although  the  noise  level  a  must  be  known  or  estimated 
in  this  improved  maximum  likelihood  estimation,  the  result  is  not  greatly  influenced  by 
the  use  of  the  exact  or  an  approximate  value:  both  give  the  same  leading  order  term  in 
the  expansion  of  the  bias. 

The  leading  order  term  can  be  removed,  however,  by  using  a  second  order  approxima¬ 
tion  to  the  log-likelihood  function.  Again,  while  this  approximation  requires  knowledge  of 
the  noise  level,  using  the  true  value  or  an  estimate  makes  little  difference  to  the  resulting 
estimate  of  the  radius. 

It  is  important  to  include  the  term  involving  curve  length  in  the  second  order  estimate, 
however:  some  simple  calculations  show  that  if  this  is  dropped  the  resulting  estimate  is 
then  even  worse  than  that  produced  by  straight  orthogonal  distance  regression.  This  is 
not  unexpected:  as  noted  in  section  4  the  linear  correction  term  tends  to  favour  points 
lying  on  the  outside  of  a  convex  curve.  Given  that  the  stochastic  process  underlying  the 
formation  of  observations  will  produce  more  points  outside  the  circle  than  inside,  in  the 
absence  of  a  term  favouring  shorter  curves  the  linear  correction  term  will  simply  reinforce 
the  tendency  of  straight  orthogonal  distance  regression  to  overestimate  the  true  radius. 

Finally,  the  results  in  this  section  show  that  the  tendency  of  orthogonal  distance  re¬ 
gression  to  underfit  noses  is  not  due  to  the  convexity  of  the  nose  per  se.  Rather  the  change 
in  curvature  round  the  nose  is  the  key  factor  that  ensures  that  more  points  are  likely  to 
be  observed  inside  the  nose  than  outside.  In  a  nose  of  constant  curvature  (i.e.,  a  circle), 
observations  are  more  likely  to  fall  outside  the  “corner”  than  inside. 


6  Examples 

In  this  section,  we  present  the  results  of  a  number  of  simulations  of  fitting  ellipses  to 
noisy  data.  The  data  points  were  generated  by  adding  identically  distributed  Gaussian 
noise  to  the  x  and  y  components  of  points  that  were  uniformly  sampled  along  the  arc 
length  of  a  given  ellipse.  To  keep  matters  simple,  the  ellipse  was  centered  on  the  origin 
and  aligned  with  the  x-y  axes,  so  that  the  curve  c  was  given  by 

y^ 

c  =  {f(x,y,a,b)  =  1  :  f(x,  y;  a,  b)  =  +  jp }■ 
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The  noisy  points  from  the  ellipse  were  fitted  by  both  orthogonal  distance  regression 
and  minimising  the  log-likelihood  function  (13)  under  the  assumptions  that  the  variance  is 
known.  In  each  of  the  two  test  cases  presented  here,  5000  points  were  sampled  uniformly 
along  the  arc  length  of  a  specified  ellipse.  Identically  and  independently  distributed  zero- 
mean  Gaussian  noise  with  a  =  0.05  was  added  to  each  of  the  x  and  y  coordinates  of  the 
sampled  points.  In  the  first  test  case,  the  parameters  of  the  ellipse  were  a  =  1  and  b  —  0.2 
(figure  3(a)),  and  in  the  more  difficult  second  test  case  they  were  a  —  1  and  b  =  0.1 
(figure  4(a)). 


The  distances  of  the  noisy  points  from  the  curve  required  to  compute  the  sum  of 
squared  orthogonal  distances  (2)  and  the  log-likelihood  function  (13)  for  particular  ellipse 
parameter  values  a  and  b  was  found  by  computing  the  closest  point  on  the  ellipse  using 
one  of  the  techniques  in  [17].  The  arc  length  and  radius  of  curvature  of  a  point  on  an 
ellipse  is  well  known  and  can  be  obtained  from  many  standard  reference  texts. 


The  contours  of  the  sum  of  squared  orthogonal  distances  and  second-order  approxima¬ 
tion  to  the  log-likelihood  function  as  they  vary  with  a  and  b  are  shown  in  figures  3(b)  and 
3(c)  respectively  for  the  first  example,  and  figures  4(b)  and  4(c)  respectively  for  the  sec¬ 
ond  example.  The  bias  in  the  orthogonal  distance  regression  estimates  is  clearly  evident. 
In  fact,  in  the  more  difficult  case  (figure  4(b)),  the  minimum  of  the  orthogonal  distance 
regression  metric  does  not  have  a  well  defined  minimum.  In  contrast,  the  minimum  of  the 
second-order  approximation  to  the  log-likelihood  function  (13)  is  better  defined  and  has 
a  bias  that  is  an  order  of  magnitude  smaller,  especially  in  the  more  difficult  example  of 
a  =  1  and  b  =  0.1. 


As  argued  earlier,  orthogonal  distance  regression  should  overfit  convex  curves,  and 
indeed  the  bias  exhibited  is  towards  uniformly  larger  ellipses.  Unfortunately  this  bias 
is  most  pronounced  at  the  ends  with  maximum  curvature,  contradicting  the  previous 
argument  for  underfitting  of  noses.  This  shows  the  strengths  and  weaknesses  of  intuitive 
arguments.  While  both  are  plausible  they  actually  argue  for  effects  in  opposing  directions, 
and  in  the  present  case  it  seems  likely  the  bias  is  actually  due  to  other  aspects  of  the 
ellipses’  geometry.  Nevertheless  bias  clearly  exists;  equally  clearly  it  is  greatly  reduced 
under  the  analysis  and  algorithm  proposed  here. 


Finally  we  make  one  further  observation  on  curve  fitting  based  on  our  experience  in 
constructing  these  examples.  The  large  sample  sizes  and  noise  levels  in  the  examples 
severely  tested  all  the  algorithms  used  for  finding  closest  points  on  the  curve:  their  odd 
failure  for  particular  configurations  in  turn  induced  large  errors  in  the  calculated  fit.  For 
example,  the  closest  point  can  be  characterised  as  the  root  of  a  quartic,  but  care  must  be 
taken  in  constructing  this  quartic  otherwise  it  becomes  degenerate  or  close  to  degenerate 
in  certain  configurations.  This  in  turn  means  the  wrong  root  is  chosen  as  the  closest 
point,  producing  an  outlier  in  the  distribution  of  orthogonal  distances.  As  is  usual  in  least 
squares  regression,  such  outliers  can  strongly  influence  the  computed  fit:  when  it  occured 
this  effect  overrode  the  bias. 
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7  Conclusions 

The  analysis  here  confirms  intuition  and  experience  by  showing  that  orthogonal  dis¬ 
tance  does  indeed  give  biased  fits  when  the  noise  level  is  significantly  large  compared  to 
the  radius  of  curvature.  Moreover,  while  maximum  likelihood  estimates  based  on  non¬ 
linear  models  are  generally  biased  [2],  the  bias  here  is  qualitatively  different  in  that  it  is 
independent  of  the  number  of  data  points.  This  makes  orthogonal  distance  regression  an 
inconsistent  estimator. 

This  has  been  noted  before  [6],  but  the  analysis  proposed  here  to  account  for  it  appears 
to  be  novel.  In  [6]  the  curve  is  defined  implicitly  and  bias  is  reduced  by  making  a  second 
order  correction  to  the  constraint  rather  than  the  penalty  function.  Minimal  justification 
is  given  for  this  step,  and  the  resulting  expression  is  not  easily  expressed  in  the  terminology 
used  here  (in  practice  it  appears  to  result  in  a  minimisation  similar  to  that  in  (15)  but 
without  the  term  involving  the  length  of  the  curve). 

The  major  advantage  of  the  analysis  here  is  that  it  highlights  the  likely  cause  of  bias 
and  inconsistency.  The  maximum  likelihood  estimation  problem  (1)  defined  by  the  stan¬ 
dard  model  does  not  have  a  fixed  set  of  parameters:  instead  the  number  of  unknown 
parameters  yn  increases  proportionally  with  the  data.  Thus  standard  asymptotic  results 
on  convergence  of  the  estimator  do  not  apply.  By  extending  the  basic  model  we  are  able 
to  integrate  out  the  unknown  yn  at  the  cost  of  an  extra  assumption:  that  the  samples  are 
uniformly  distributed  along  the  true  curve.  The  result  is  a  more  complete  description  of 
the  problem  that  includes  all  its  stochastic  components  and  requires  the  modeller  to  make 
explicit  assumptions  about  each.  This  focus  on  a  statistical  model  for  curve  fitting  has  the 
added  benefit  in  that  it  provides  a  framework  for  subsequent  processing:  for  example,  the 
problem  of  merging  separate  curve  segments  can  be  veiwed  as  an  exercise  in  hypothesis 
testing  (the  tests  may  be  based  on  comparisons  of  the  distributions  of  the  residuals  from 
separate  fits  vis  a  vis  the  residuals  of  a  joint  fit). 

There  are  a  number  of  obvious  issues  and  extensions  that  we  do  not  have  the  space 
to  address  in  detail  here,  but  that  deserve  mention.  First,  in  many  problems  the  curve 
length  |c|  cannot  be  expressed  in  closed  form  but  must  be  approximated.  Ideally  this 
would  be  computed  directly  by  some  form  of  interpolation  of  the  calculated  closest  points 
x*.  Second,  in  practice  the  distribution  of  samples  along  the  true  curve  is  not  likely  to 
be  uniform,  especially  if  the  curve  is  infinite.  In  this  case  the  model  could  be  extended 
to  include  the  unknown  density  and  this  could  be  estimated  as  part  of  the  likelihood 
maximisation.  (This  would  also  remove  the  need  to  estimate  |c[.)  Third,  the  data  points 
and  their  perturbations  from  the  true  curve  are  often  higly  correlated:  this  occurs  whenever 
a  curve  is  fitted  to  an  extracted  edge  or  skeleton.  In  this  case  the  model  should  be  extended 
to  include  a  model  for  correlated  noise  and  estimation  of  any  undetermined  parameters  in 
this:  if  it  isn’t,  in  our  experience  the  fits  produced  are  still  reasonable,  but  post  processing 
activities  such  as  curve  merging  [16]  fail  because  the  statistical  assumption  on  which  they 
are  based  are  violated. 

Finally  the  most  important  question  is  in  what  circumstances  is  the  fit  obtained  by 
minimising  (13)  likely  to  produce  better  results  than  standard  orthogonal  distance  regres¬ 
sion?  While  we  do  not  have  a  definitive  answer,  three  observations  on  it  are  worth  noting. 
First,  in  practice  there  no  significant  difference  in  the  computational  complexity  of  the 
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two  minimisations.  The  computational  cost 

point  °n  the  “"w  n^ic  pZos^  Le  costs  on  the 

°f  ‘ ^nUonThatTho  underlying  true  samples  are  distributed  uniformly  along  the  curve, 
assumption  that  tne  unuciry  h  information  to 

The  assumption is e „  likely  to  be  same 

the  contrary,  but  if  it  is  not  true  the  «„ffeests  that  bias  will  dominate 


radius  of  curvature. 
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Figure  Al:  A  typical  ionogram.  Frequency  is  plotted  along  the  x  axis  and  signal  travel-time 
along  the  y-axis.  The  grey-level  indicates  the  strength  of  the  received  signal. 

Appendix  A:  Autoscaling  Ionograms 


High  frequency  (HF)  radio  communications  over  long  distances  rely  on  reflections  from  the 
ionosphere  to  bounce  the  wave  from  transmitter  to  receiver  around  the  curve  of  the  earth. 
The  capacity  of  this  indirect  channel  is  a  (complicated)  function  of  transmission  frequency 
and  the  state  of  the  ionosphere,  which  changes  continuously.  Therefore  communications 
systems  and  over-the-horizon  (OTH)  radar  systems  that  make  use  of  ionospheric  reflection 
require  up-to-date  information  on  the  state  of  the  ionosphere.  This  information  can  be 
gleaned  from  ionograms:  these  are  visual  displays  of  suitably  chosen  test  signals.  Figure  Al 
shows  a  typical  ionogram:  the  horizontal  axis  indicates  the  frequency  of  the  transmitted 
signal,  the  vertical  axis  indicates  the  time  taken  for  the  signal  to  reach  the  receiver  and 
the  grey-level  indicates  the  intensity  of  the  received  signal. 

For  the  information  contained  in  figure  Al  to  be  useful  for  HF  communications  and 
OTH  radar,  is  has  to  be  scaled,  i.e.,  it  has  to  be  interpreted  by  a  trained  ionospheric 
physicist.  This  involves  identifying  the  various  curves  or  traces  in  the  image  along  with 
their  features,  and  associating  these  with  paths  taken  by  the  signal  and  the  atmospheric 
structure  on  those  paths  (refraction  by  the  ionosphere  varies  with  frequency  plus  a  signal 
path  may  include  more  than  one  reflection  between  transmitter  and  receiver).  This  infor¬ 
mation  can  be  used  both  to  select  optimal  operating  frequenices  for  communications  or 
radar  systems  and  to  infer  the  structure  of  the  ionosphere.  Ionospheric  sounding  networks 
are  now  recording  ionograms  over  a  multitude  of  paths  at  intervals  of  minutes,  however, 
and  the  resulting  data  set  far  outstrips  the  supply  of  physicists  available  to  interpret  it. 
Hence  there  is  increasing  need  for  an  autoscaling  system  that  can  automatically  interpret 
recorded  ionograms. 

Development  of  a  prototype  autoscaling  system  [14,  15,  16,  17]  threw  up  a  number  of 
significant  curve  fitting  and  post-processing  problems  which  in  turn  motivated  the  work 
presented  here.  An  autoscaling  system  must  be  able  to  identify  traces  and  fit  them  by 
a  simple  parametric  model  that  can  then  be  used  for  purposes  such  as  classification  or 
tracking  of  traces.  While  orthogonal  distance  regression  is  an  obvious  and  useful  starting 
point  for  such  a  fitting  procedure,  various  features  of  the  autoscaling  problem  require  a 
more  sophisticated  model  of  how  data  is  generated  and,  consequently,  of  the  optimal  fitting 
process. 

For  example,  any  trace  fitting  routine  must  give  an  accurate  estimate  of  the  nose  of 
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,,  ,nTnin_nt  trace  as  this  gives  the  optimal  frequency  for  transmission  (transmissions  at 
lower  frequencies  will  be  corrupted  by  multi-pathing;  transmissions  at  higher  frequencies 

l  ..  '}  t  d  bv  the  ionosphere).  Likewise  the  particular  measurements  selected 

for  uso  i,>  trace  BHi„g  arc  often  correlated  or  have  other  significant  statistical  structure 
V  Tv  fit  ed  traces  must  often  be  further  processed,  c.g„  separate  segments  of  the  same 
‘  nmst  be  merged.  All  of  these  activities  require  a  good  statistical  model  of  trace 
r  f  that  can  then  be  used  to  construct  well-founded  fitting  and  merging  routines. 
T“cr  aims  to  provide  a  statistical  setting  for  orthogonal  distance  regression  that 
enables  it  to  be  extended  to  meet  these  demands. 
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Appendix  B:  Moments  of  the  Conditional 
Distributions  Associated  with  Circles 


This  appendix  derives  various  closed  form  expressions  and  asymptotic  expansions  for  the 
moments  of  the  conditional  distributions  of  noisy  observations  of  points  uniformly  dis¬ 
tributed  along  circles. 

Let  c  be  a  circle  of  radius  r  centred  on  the  origin.  We  wish  to  calculate  the  moments 
of  the  conditional  distribution  p(x|c,  a)  of  (16),  in  particular  the  moments  E[pk ]  where 
p  =  j | x 1 1 .  From  (16), 

E[pk]  =  ±  J™  pke-^l2a\-rpl^h(^)pdp 

■  ^r^*(SS)* 

=  (2a2)‘/V^r(l  +  l)  m  (1  +  1.1,^)  (Bl) 

from  6.631.1  of  [7].  Trivially,  E[p°]  =  1.  To  compute  E\p],  note  that  9.215.3  and  9.212  of 
[7]  give 

l.  jFi  ( 1  /  2 , 1,  z)  =  ez/2I0{z/2), 

^(1/2,1  ,z)  =  x  ii?i(3/2, 2,  z), 
az  2 

« iFi(3/2, 2, z)  =  iFi(3/2,l,z)  -  iFi(1/2,1,^). 

Therefore, 

^  1  ^ 

iFi(3/2,  l,z)  =  2*  —  iFi(-,M)  +  iFit^M) 

=  +  ^  + 


Substituting  this  back  into  (Bl)  gives 


£[p]  = 


7T  r2/4ff2 
2 


+  l  /o 


+  2^Jl 


We  are  especially  interested  in  the  asymptotic  region  in  which  r/cr  »  1:  in  this  case 
substituting  in  the  appropriate  asymptotic  expansions  for  I0  and  I\  gives 


A  closed  form  for  E[p2}  follows  from  noting  that  9.215.1  and  9.212.4  of  [7]  gives 
iFi(l,M)  = 

iFi(2,  1,  z)  =  (z  +  l)iFi{l,l,z)^{z  +  l)ez. 
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Substituting  this  back  into  (Bl)  gives 

E[p 2]  =  r 2  +  2ct2. 

Next  E[p4}  follows  from  noting  that  applying  9.212.4  of  [7]  again  gives 

(z2  +  4z  +  2)  , 
iFi(3, 1,  z)  =  - x - e 


so  that 

E[p 4]  =  r4  +  8r2cr2  +  8(t4. 


(B3) 


(B4) 
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Appendix  C:  Fitting  a  Surface  to  Points  in 

Space 


This  section  extends  the  results  in  the  body  of  the  paper  on  fitting  a  curve  through  points 
in  a  plane  to  the  problem  of  fitting  a  d- 1  dimensional  surface  S  through  a  set  of  points 
{xn}£Lx  in  .  In  particular  it  derives  a  second  order  approximation  to  the  log-likelihood 
function  based  on  a  second  order  asymptotic  approximation  to  the  integrals  defining  the 
conditional  probabilities  p(x„|iS',  a).  It  is  not  as  easy  to  work  with  the  quadratic  surfaces 
that  are  natural  local  second  order  approximations  to  an  arbitrary  surface  as  it  was  with 
circles.  Therefore  we  instead  develop  a  local  approximation  to  S  about  the  closest  point 
x*  to  xn  as  a  parabolic  surface  whose  form  and  orientation  is  determined  by  the  principal 
radii  of  curvature  of  S  at  x* .  The  approach  can  be  readily  extended  to  the  more  general 
problem  of  fitting  an  d-k  dimensional  surface  5  through  points  in  Rd ,  but  for  simplicity 
we  shall  limit  ourselves  to  the  case  k  =  1. 

In  particular,  let  S  be  a  given  d—  1  dimensional  surface  in  and  let  x  be  a  noisy 
observation  of  S.  As  before  we  assume  that  the  probability  that  x  is  a  perturbation  of  a 
point  y  on  S  is  given  by 

p(x|y’Si")=id^e"l|x"yl,V2"Sdx' 

Thus  the  a  priori  assumption  that  y  will  be  uniformly  distributed  over  S  gives 


p(x|£>) 


[2ira2]d/2  |S|  Jy{s)eS 


I 

J  vf 


-||x-y(S)||2/2<72ds 


(Cl) 


where  |5|  is  the  area  of  S,  and  ds  is  an  element  of  surface  area. 

We  now  set  up  a  local  approximation  to  the  integral  in  (Cl).  To  do  so,  we  first  define 
a  local  coordinate  system  as  follows.  Let  x*  be  the  closest  point  to  x„  on  S:  x*  will  be 
the  origin  for  the  coordinate  system.  Next  let  Tn  be  the  tangent  plane  to  S  at  x*.  The 
directions  associated  with  the  principal  curvatures  of  S  at  x*  form  an  orthogonal  basis 
for  Tn  [9];  let  the  unit  vectors  in  these  directions  be  ei, . . .  ,  e^-i-  Then,  since  xn  —  x*  is 
orthogonal  to  Tn,  we  may  take  the  final  coordinate  direction  ed  to  be  the  unit  vector  in 
the  direction  xn  —  x* . 

As  before,  let  dn  =  ||xn  —  x*  ||.  Rather  than  interpreting  dn  as  a  signed  distance  and 
using  this  to  describe  the  local  geometry,  however,  we  now  instead  take  the  signs  of  the  cur¬ 
vatures  ,  Kd- 1  to  be  with  respect  to  the  direction  ed-  Thus  in  this  coordinate  system 

the  point  xn  has  coordinates  (0, . . .  ,  0,  <in),  and  the  surface  S  can  be  locally  represented 
as  the  following  parabola  over  the  tangent  plane 

S  ~  {(?/i,  -  -  -  ,yd-i,yd(yi,---  ,yd-i))  ■  (yi,--- ,yd-i)  £Tn} 


where 


DSTO-RR-0242 


Likewise  we  can  approximate  the  surface  element  ds  by  the  elemental  product  dy\ . . .  dyd-i- 
Thus  (Cl)  can  be  approximated  by 


p(x.n\S,o) 


[27 


,21^/2 


\S\ 


j 

JBd-l 


-\\xn-y(s)\\V2^dyi  'dyd_ 


Next  expanding  the  distance  ||xn  -  y(s)||2/2cr2  as 

fd-i 

2a2"  ~n  JK"'"  2  a2 


1  '|x„  -  y(.s)||2  =  WZ2  I  S  Vi  +  “  d») 


1 

2a2 

1 
2 


i— 1 
rf-1 


X>*2+ 

t=l  V  i=l  / 


and  making  the  change  of  variables  z\  =  yj/a  gives 
-dl/2<T2 


p(xn|5,  cr) 


/  exp 

-l(E-rf)2 

•  exp 

r  i  d~i  i 

-  “  Kidn)zi 

/Rd-1 

.  Z  i=l 

dzi 


[27r]d/2a|5| 

Under  the  assumptions  that 

dji 

ok;  <C  1,  dnKi  <C  1,  and - 1 

o 

then  an  asymptotic  expansion  of  the  integral  above  in  each  variable  gives  to  first,  order 

P-Snj2a2 

p(xn|5,  o)  ~  m  -  Mn)_1  2- 

v2tt  |5|(7  fJi 


Thus  the  likelihood  function  is 


Z(x„|5,a)  =  -lnp(x„|5,<r 
~  hi 


(C2) 


In 


In 


\/27r  1 5 1 a  +  I^E 1,1(1  “M"> 

i= 1 

*2  A  ^ 

’’  "E-- 


\/27r  (S|cr 
\piA  |5|a 


+ 


2a2 


2  —  1 


+ 


d2  (d  -  1  )dr 


2  a2 


where  Ki  is  ^ie  mcan  curvature  of  5  at  x*.  Finally,  given  N  independent 

observations  the  log-likelihood  function  is  approximated  by 


N 

Z(x i, . ..  ,xN:S,o)~J2 

n— 1 


{d  l)tindn 

+  JV 

In  |S|  +  In  \Z2tto 

_2  a2  2 

L  J 

•  dzd-i- 
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This  approximation  is  an  obvious  extension  of  (13)  after  making  the  simple  replacement 
of  the  curvatures  m  by  the  corresponding  radii  of  curvature  r*  =  l//q.  Again,  if  the  vari¬ 
ance  is  unknown  it  can  be  eliminated  by  observing  that  at  the  minimum  o  jv  Sn=i  • 
Substituting  in  this  value,  exponentiating  and  squaring  gives  that  an  approximation  to 
the  maximum  likelihood  estimate  of  S  will  be  the  minimiser  of 

(jv  \  "  d  i  n 

4 £ di )  ■  exp  £ Kndn 

n=l  /  L  n= 1 
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